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EFFICIENT COMPUTATION OF THE BERGSMA-DASSIOS SIGN 

COVARIANCE 

LUCA WEIHS, MATHIAS DRTON, AND DENNIS LEUNG 


Abstract. In an extension of Kendall’s r, Bergsma and Dassios (2014) intro¬ 
duced a covariance measure r* for two ordinal random variables that vanishes 
if and only if the two variables are independent. For a sample of size n, a di¬ 
rect computation of t* , the empirical version of t* , requires O(n^) operations. 
We derive an algorithm that computes the statistic using only 0(n^ log(n)) 
operations. 


1. Introduction 

Kendall’s r (Kendall, 1938) and Spearman’s p (Spearman, 1904) are popular 
measures of dependence between two random variables X and Y. However, both 
have the undesirable property that they may be equal to zero even when X and Y 
are not independent. Addressing this weakness, Bergsma and Dassios (2014) have 
defined a new coefficient, t*, which, under mild conditions on the joint distribution 
of {X, K), is zero if and only if X and Y are independent. However, a computational 
price is to be paid for this property as a naive computation of t*, the empirical 
version of t* , requires O(n^) time for a sample of size n. 

In this paper we present an algorithm which computes t* in 0{n^ log(n)) time, in¬ 
spired by a similar improvement for computing (the empirical version of) Kendall’s 
T. Indeed, by leveraging binary tree algorithms and observing that Kendall’s statis¬ 
tic depends only on the relative order of data points, Christensen (2005) showed 
that Kendall’s r can be computed in 0(n log(n)) time rather than O(n^). We follow 
a similar strategy by exploiting the fact that computing t* relies only on the relative 
ordering of quadruples of points. Due to excessive time requirements, Bergsma and 
Dassios limit their computational examples to sample sizes with n < 50 and suggest 
approximating t* by random subsampling for larger samples. As will be shown in 
Section 4, our algorithm computes t* exactly in less than a second for sample sizes 
in the thousands. 


1.1. Background and Setup. Given a sample {xi,yi),{xn,yn) of points in 
define the statistic 


( 1 . 1 ) 


t* := 


(n — 4)! 
n! 


X! a{xi,Xj,Xk,Xi)a{y^,yj,yk,yi) 

distinct 


where 


a{zi,Z2,Z3,Z4) ■■= sign(| 2 ;i - Z 2 I -I- |z 3 - Z 4 I - jzi - Z 3 I - |z 2 - Z 4 I). 


Date: January 2, 2016. Key words and phrases. Binary tree, Kendall’s tau, nonparametric 
correlation, Spearman’s rho, rank correlation, test of independence. 
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(a) Concordant (b) Discordant 

Figure 1. Relative position of points within quadrants does not 
matter, only that they remain in their respective quadrants. 


Here t* is the U-statistic, U standing for unbiased, corresponding to the popula¬ 
tion coefficient t* := Ea{Xi, X 2 , X 3 , X 4 )a(Yi,Y 2 ,Y 3 ,Y 4 ) of Bergsma and Dassios 
(2014) where (Xi, Yi),..., (X 4 , 14 ) are random vectors drawn independently from 
some bivariate distribution on Bergsma and Dassios (2014) introduce t* not as 
a U-statistic but as the closely related biased V-statistic; we consider the U-statistic 
as it simplifies some of the computations in Sections 2 and 3 but present modifi¬ 
cations to our algorithm that allow one to compute the V-statistic in Appendix A. 
A comprehensive overview of U/V-statistics and their properties can be found in 
Serfling (1980). 

As noted by Bergsma and Dassios (2014), we may rewrite the function a as 
a{zi,Z2,Z3,Z4) = I{zi,Z3 < Z 2 , Z 4 ) +I{zi,Z3 > Z 2 , Z 4 ) 

- I{zi,Z2 < Z3,Z4) - I{zi,Z2 > Z3,Z4) 

where I(zi,Z2 < 23 , 2 : 4 ) is the indicator of the event max(zi,Z2) < niin( 2 : 3 , Z4). 
After rewriting a in this way we see that computation of the t* statistic requires 
only knowledge of the relative positioning of the observations for which we make 
the following definitions. Let (xi,yi), ..., (3:4,1/4) be four points relabelled so that 
xi < X2 < X3 < X4. We then say that the points are 

,, .p X2 = X3 or there exists a permutation tt of {1, 2, 3,4} 
mseparabte it ,, , ^ 

so that 2 /^( 1 ) < ?/^( 2 ) = 2 /^( 3 ) < 21 , 7 ( 4 ), 

and if they are not inseparable, then we call them 

f concordant if max(2/i, 2 / 2 ) < niin(2/3, 2 / 4 ) or max(2/3, 2 / 4 ) < niin(2/i, 2 / 2 ), 

} discordant if max( 2 /i, 2 / 2 ) > min(2/3, 2 / 4 ) and max(2/3,2/4) > inin(2/i, 2 / 2 )- 

These definitions categorize all quadruples of points, that is, any quadruple of points 
must be exactly one of inseparable, concordant, or discordant. Moreover, when all 
coordinates are distinct any collection of four points will be either concordant or 
discordant, see Figure 1. We motivate calling points inseparable by noting that, in 
the X 2 = X 3 case, we cannot draw a line parallel to the 2 /-axis that separates the 
X values into two groups. Similarly in the case of 2/77(1) < 2/77(2) = 2/77(3) ^ 2/77(4), 
there exists no such line parallel to the x-axis that separates the y values into two 
groups. 
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We will derive two algorithms for the computation of t*, the first works only in the 
case that the data contains no ties, that is all Xi,x„ are distinct and similarly for 
yi,and the second works for all data. While the second algorithm is strictly 
more general than the first it is also substantially complicated by the need to 
consider the case of inseparable points. We present the algorithm for data without 
ties in Section 2 and give the general algorithm in Section 3. 


1.2. A Preliminary Lemma. Before moving on, it will be useful to rewrite t* 
to capture a certain permutation invariance and state a basic, but very useful, 
lemma. Let C(n, 4) = : l<i<j<k<l< n}, and S 4 be the 

set of permutations on 4 elements. For ease of notation, for any tt € S 4 and 
(zi, Z 2 , Z 3 , Z 4 ) S K'* we define Z 7 r(i, 2 , 3 , 4 ) := (-Z 7 r(i): ■■■,' 2 ^ 7 r( 4 ))- We may then rewrite 
( 1 . 1 ) as 


t* = 


(1.3) 


where 


(n — 4)! 


(n — 4)! 


n! 


E E 

{2,_7,/c,/}GC(n,4) 7t^S4 

E/ bijki, 


^ijkl ■ ^ ^ ) 

is clearly invariant to any permutation of A:, L 

We now characterize the possible values bijki may take. 

Lemma 1 . 1 . Lef A = {(xi, yi), (x2, y2), (2^3,2/3), ( 2 ^ 4 , J/ 4 )} C Then 

{ 16 if the points in A are concordant 

—8 if the points in A are discordant 

0 if the points in A are inseparable 

The proof of Lemma 1.1 is a straightforward but lengthy case-by-case analysis 
and we defer it to Appendix B. 


2 . The Algorithm for Data Without Ties 


Throughout this section we assume that (xi, yi),..., (x„, y„) contain no ties, that 
is, xi,...,x„ are pairwise distinct and so are yi,...,y„. As there are no ties, every 
quadruple of points is either concordant or discordant. It follows from Equation 
(1.3) and Lemma 1.1 that 


t* = 


(n — 4)! 


E 


ijkl 


{i,j,k,l}^C {71,4) 


(n — 4 )! 


E 16 • ..., ixi,yi) are concordant}) 


{i,j,k,l}GC{n, 4 ) 


- 8 • /({(xi,y 4 ,..., (x/,yi) are discordant}) 


(n - 4V 

= ^-E(16-A,-8-Ad) 

71.1 


(n — 4)! 
nl 


(24 • A,) 


1 

3 ’ 


( 2 . 1 ) 
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where Nc and Nd are the numbers of concordant and discordant quadruples in 
ixi,yi), ixn,yn), respectively, and the last equality holds since every quadruple 
of points is either concordant or discordant implying that (") = Nd + Nc- Thus 
computing t* requires only computing the number of concordant quadruples of 
points. We now show that this can be done efficiently. 

Suppose we have relabeled the points so that xi < X 2 < ■■■ < Xn- Rewriting 
sums we have that 

Nc= ^ I{{xi,yi), ixj,yj), {xk,yk), {xuyi) are concordant)) 

l<i<j <k<l<n 

= E E E I{{x^,yi),{xj,yj),{xk,yk)Axhyi) are concordant)) 

3<fc<n—1 k<il<n <k 

= E E E Hvi^yj < yk,yi) + iiyk,yi < 

3<fc<n—1 k<il<n <k 



where we define 

M^{k,l) := |{z : 1 < i < k, yi < uim{yk,yi)}\, 

M^{k,l) := |{i : 1 < i < k, yi > max(j/fc,y;)}|. 

The last line in the above summation is, effectively, the algorithm. Note that the 
summation is over 0{n?) terms and, consequently, if we can find M^(k,l) and 
My{k,l) in 0(log(n)) time then we have found an algorithm for computing 
in 0(n^ log(n)) time. To find M^{k,l) and My{k,l) in 0(log(n)) time we use 
a binary tree data structure with an appropriate balancing algorithm to ensure 
that inserts and searching can be done in 0(log(n)) time. One example of this 
type of data structure are the so-called red-black trees (Guibas and Sedgewick, 
1978). In particular, given that we have inserted the values yi,y 2 , ■■■,yk-i into 
a red-black tree we may insert another yk into the tree in 0(log(fc)) time and a 
simple extension of the traditional red-black framework allows one, for any y, to 
find \{l < i < k — 1 : yi < y}\ and \{1 < i < k — 1 ■. yi > y}\ \n 0(log(fc)) time. 

Combining the above observations. Algorithm 1 gives an 0(n^ log(u)) procedure 
for finding the number of concordant quadruples which is easily extended to a com¬ 
putation of t* via Equation (2.1). Note that in Algorithm 1 there is a preprocessing 
step in which we sort the xi,...,Xn values in ascending order and then reorder the 
yi to match this new order. Since this preprocessing can be done in worst case 
0{n\og{n)) time with a number of algorithms, merge-sort for example, it is not a 
significant component of the overall asymptotic run time analysis. 

3. The General Algorithm 

Now suppose that there are no restrictions on the values of {xi,yi ),..., {xn,yn) 
and that we have reordered the points so that xi < ... < Xn- For any 3 < k < I < n, 
let 

(3.1) top{kJ) = \{1 < i < k : x^ ^ Xk and yt > max(yfe, ?//)}|, 

(3.2) mid{kj) = \{1 < i < k : x^ ^ Xk andmin(yfc,y;) < yi < max(yfc,y/)}|, 

(3.3) bot(k,l) = \{l < i < k -. Xi ^ Xk and yi < min(yfe, y;)}|. 
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Algorithm 1 

1 : procedure NuMCONCORDANT((a;i, ?/„)) 

2 : X ^ {xi,....,Xn) 

3: y ^ {yi,--,yn) 

4: Sort X in ascending order and relabel y to match this new order 

5: rbTree ^ empty red-black tree 

6: totalConcordant ■<— 0 

7: for k = 1,n — 1 do 

8: for £ = k + 1,n do 

9: minY ^ min{yk,ye) 

10 : maxY <—max{yk,ye) 

11: numLess number of elements < minY in rbTree 

12: numGreater ^ number of elements > maxY in rbTree 

13: totalC oncordant = totalC oncordant + + ^numGpater^ 

14: Insert y/^ into rbTree 

15: return totalC oncordant 


eqMax 


eqMin 


top 

1 

1 

1 

1 




* 

mid 

1 

1 

1 


1 

bot 

1 

1 

1 

1 


Figure 2 . Partitioning of the points with x value strictly less than 
two given points. Solid lines correspond to eqMax and eqMin, the 
points whose y-values equal the maximum or minimum of the y 
values of the two given points, respectively. 


(3.4) eqMin{k,l) = \{1 < i < k : x^ ^ Xk and yi = inm{yk,yi)}\, 

(3.5) eqMax(k,l) = \{1 < i < k ■. Xi ^ Xk and yi = max( 2 /j,, ?//)}|. 

These quantities correspond to a partitioning of the points {xi,yi) with i < k and 
Xi ^ Xk- We illustrate this partitioning in Figure 2. For fixed 3 < k < I < n we 
have, by Lemma 1.1 and since xi < ... < x„, 

= 16 • |{1 < i < j < fc : i,j, fc, I correspond to concordant points}| 

— 8 • |{1 < i j < fc : i,j,k,l correspond to discordant points}] 

'-V-' 

:=Adidfc.O 

= 16-lVcon(fc,0-8-^dis(fc,0- 
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Hence, similarly as in the case without ties, we may write 

= y]] y]] y]] ^ijM 

{ij,/c,Z}GC(n,4) 3<fc<n—1 k<l<n l<i<j<k 

= 16-iVeo„(fc,0-8-A^dis(fc,0- 

3<fc<n—1 k<l<n 


Again the last line of the above summation is effectively the algorithm. Since 
the sums are over O(n^) terms, if we can show that Ncon{k,l) and A^dis(fc, 0 can 
be computed in 0(log(n)) time then we have obtained an 0(n^log(n)) algorithm 
for computing t*. We show next that this is indeed possible, beginning with the 
observation that 


(3.6) 

if Vk = yi then 

(3.7) 


Ncon{kJ) 


^op{k, Z)^ ^ot{k, 0^ 


Ndis{kj) = 0, 


and if yk ^ yi then 

(3.8) 


Ndis{k, 1) = topik, 1) ■ {mid{k, 1) + eqMin{k, 1) + bot{k, 1)) 
+ bot{k, 1) ■ {mid{k, 1) + eqMax{k, 1)) 

+ eqMin{k, 1) ■ {nnid{k, 1) + eqMax{k, 1)) 

+ eqMax{k, 1) ■ mid{k, 1) 


mid{k, 1) 
2 


E 

yGunique{k,l) 


\{l < i < k : Xk ^ Xi and yi = y}\ 
2 


where 


unique{k,l) := {y^ : 1 <i<k and Xt ^ Xk and mm{yk,yl)<y^<max{yk,yl)}■ 


Suppose we have a red-black tree into which we have inserted all yi with 1 < z < fc 
and Xi ^ Xk- Then it is clear that the quantities in Equations (3.1)-(3.5) can 
each be computed in 0(log(fc)) time. Note that, unlike in the untied case, we 
require that the red-black tree not include any yi values corresponding to Xi = Xk', 
accomplishing this algorithmically is very simple: as we iterate across the Xk values 
we delay inserting their associated yk values into the red-black tree until we reach 
a xi with xi ^ xi-i, upon reaching such an xi we insert all postponed y values into 
the red-black tree and then restart the postponing of y values starting with yi. 

We see that, as in the discussion of Algorithm 1, we can progressively compute 
almost all of the quantities in Equations (3.6) and (3.8) with each iteration taking 
0(log(n)) time. The only complication is the computation of 


(3.9) 


y^unique{k,l) 


/|{1 < i < k : Xk Xi and yi = z/}| 

V 2 


which corresponds to all quadruples of points {xi,yi), {xj,yj), {xk,yk), ixi,yi) for 
which min(z/fc, z/i) < yi = yj < max^yk^yi). These are inseparable and are being 
over-counted by Note that this summation is in the reverse order of 

what we would like in order to simply generalize Algorithm 1. In particular, there 
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is a condition on the y values corresponding to i and j which is suppressed by the 
aggregate counts available from a query on a red-black tree. We have, however, 
already established a methodology to count values such as (3.9). In particular, 
note that 


E E 

3<fe<Z<n y^unique{k,l) 


= ^ I{{xk ^ Xj and yi = yj and inin{yk,yi) < yi < max(yfc, j/;)}) 

l<i<j<k<l<n 

= E ^ ■ \{k : j < k < n and Xk ^ Xj and y^ < yk}\ 

l<i<j<n—2 


\{k : j < k < n and Xk ^ Xj and yj > yk}\ 

^ ■v' 

:= bof(j) 


(3.10) 


E E = Vj}) ■ top*{j) ■ bot*{j). 


It follows that all that is needed to compute the total contribution of the term in 
Equation (3.9) to t* is to run a modified version of Algorithm 1 across the data in 
reverse order. Our algorithm becomes the following: 

(i) Perform a first pass across the data where we ignore the effect of (3.9) and 
count all other quantities. 

(ii) Perform a second pass across the data in reverse order to compute (3.10). 

(hi) Appropriately combine the results of (i) and (ii) to obtain t*. 

This amounts to over-counting discordant quadruples on a first pass and then un¬ 
doing this over-counting on a second pass. Since both of these passes over the data 
require 0(n^log(n)) time, our general Algorithm 2, which leverages the above ob¬ 
servations, computes t* in 0{'n? log(n)) time. An implementation of Algorithm 2 
is available in the R package TauStar accessible via CRAN, the Comprehensive R 
Archive Network^ (R Core Team, 2015; Weihs, 2015). 


4. Simulations 

We test the run times of Algorithm 2 and a naive implementation, both written 
in C-I--I- and available through R in previously mentioned TauStar package, for 
sample sizes n ranging from 100 to 300; the implementation of Algorithm 2 uses 
the red-black tree C library of Martinian (2005). The results of these simulations 
are presented in Table 1. As the table shows, the 0(n'^) running time of the naive 
algorithm becomes already a practical concern for sample sizes in the hundreds 
while Algorithm 2 is essentially instant for such sample sizes. Table 2 provides a 
perspective on the run time of Algorithm 2 for substantially larger samples. 

It is possible to approximate C, or in other words, estimate r* by a Monte-Carlo 
subsampling procedure where, for small m < n, subsets of size m are repeatably 
selected from the data at random and the value of t* on each of these subsets is 
then averaged. Indeed, the case of m = 4 is a strategy suggested by Bergsma and 
Dassios (2014). While our algorithm makes the computation of t* on moderate 


1 


See https: //cran. r-project.org/web/packages/TauStar/index.html 





LUCA WEIHS, MATHIAS DRTON, AND DENNIS LEUNG 


Algorithm 2 

Algorithm for efficiently computing t* in the general case. Comments are displayed 
in gray._ 

1 : procedure t*{{xi,yi),...,{xn,yn)) 

2: X ^ {Xi,....,Xn) 

3: y ^ {yi,-;yn) 

4: Sort X in ascending order and relabel y to match this new order 

5: rbTree •<— empty red-black tree > Used in first pass through data. 

6 : savedYValues •<— empty list > A list to store y values whose insertion 

into the red-black tree is delayed. 

7: totalConcordant -<—0 t> Total concordant quadruples counted so 

far. 

8 : totalDiscordant -<—0 > Total discordant quadruples counted so 

far. 

9: for fc ^ 1,...., n — 1 do 

10 : if fc 7 ^ 1 and Xk-i ^ Xk then > If fc 7 ^ 1 and Xk 7 ^ Xk-i insert all de¬ 

layed y values into the red-black tree. 
In any case, save yu to be inserted in 

the tree on some future iteration. 

11 : for yV al in savedYV alues do 

12 : Insert yVal into rbTree 

13: Empty the list savedYValues 

14: Append y^ to savedYValues 

15: for I! ^ fc -I- 1,..., n do > Loop over I > k and use (3.6), 

(3.7),(3. 8 ) while ignoring contributions 
of (3.9). 

16: minY min(yf^,y() 

17: maxY ^max{yk,ye) 

18: top ■<— number of elements > maxY in rbTree 

19: mid ■‘r- number of elements < maxY and > minY in rbTree 

20 : bot ^ number of elements < minY in rbTree 

21 : eqMin ^ number of elements equal to minY in rbTree 

22 : eqMax •<— number of elements equal to maxY in rbTree 

23: totalConcordant totalConcordant + -|- 

24: if minY 7 ^ maxY then 

25: totalDiscordant totalDiscordant+[™'^^+top-mid+top- 

bot + mid ■ bot + eqMin ■ (top + mid + 
eqMax) -I- eqMax ■ (mid -I- bot) 

26: t> In the next loop we will run along the data in reverse to undo the 

over-counting resulting from ignoring the contribution of (3.9). 

27: Empty the fist savedYValues 

28: rev RbTree empty RB tree Used in second pass over the data. 

29: for j •(— n,...., 2 do 


to large samples feasible, an approximate strategy will be necessary for very large 
samples. Unfortunately, resampling procedures require choosing a number of re¬ 
sampling iterations and, as is shown by Table 3, choosing too few iterations can 
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30: if j 7 ^ n and Xj+i ^ Xj then > Inserting the delayed values similarly as 

in Line 10. 

31: for yVal in savedYValues do 

32: Insert yVal into revRbTree 

33: Empty the list savedYValues 

34: Append yj to savedYValues 

35: ior i <— j — 1, ....,1 do > Use (3.10) to compute the number of 

over counts. 

36: minY min{yi,yj) 

37: maxY ^max{yi,yj) 

38: top ■<— number of elements > maxY in revRbTree 

39: bot ■<r- number of elements < minY in revRbTree 

40: if minY = maxY then 

41: totalDiscordant •<— totalDiscordant — top ■ bot 

42: return n(^n-i){n- 2 )(n- 3 ) ’ totalConcordaut — 8 • totalDiscordant) 


Table 1. Run times of the naive algorithm and Algorithm 2 for 
various sample sizes (in seconds and averaged over 10 samples). 


Sample Size 

100 

150 

200 

250 

300 

Algorithm 2 

0.0009 

0.0023 

0.0043 

0.0072 

0.01 

Naive Algorithm 

0.287 

1.55 

5.58 

14.34 

28.95 


Table 2. Run times of Algorithm 2 for larger sample sizes (in 
seconds and averaged over 10 samples). 


Sample Size 

1000 

3250 

5500 

7750 

10000 

Algorithm 2 

0.1265 

1.7354 

5.2744 

11.1833 

19.115 


Table 3. Sample variance of resampling-based estimates of r* 
relative to the sample variance of t* computed for all data. Here 
relative variance is the ratio of the former and the latter variance. 
The variances are computed from 1000 samples of size n = 1000. 
Resampled subsets were of size m G {4,30}. The samples were 
drawn as pairs of independent N{0, 1) random variables. 


# Resamples 
Relative Var. 

(m = 4) 
Relative Var. 
(to = 30) 


200 400 800 1600 3200 6400 12800 

3932.84 2118.1 911.67 472.67 230.82 115.57 57.03 

8.24 4.19 2.43 1.67 1.21 1.11 1.04 


result in a estimator with large variance. Table 3 also suggests that a choice of 
TO > 4 may be useful.^ 

^ R code to reproduce the results of Tables 1-3 can be found on the first author’s webpage: 
http://www.stat.Washington.edu/~lucaw/public_resources/eff_coinp_2015/tables .R 
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5. Conclusion 

We have presented an algorithm which computes the U -statistic t* corresponding 
to the T* sign covariance of Bergsma and Dassios (2014) in 0(n^ log(n)) time, 
substantially outperforming a naive implementation. The computational savings in 
our algorithm are driven by the use of binary trees and the permutation invariance 
inherent in t* (recall Lemma 1.1). 


Appendix A. Modifications for the V-Statistic 

This section provides an overview of necessary modifications to Algorithm 2 in 
order to compute the V-statistic version of t*. Suppose, as usual, that we have 
reordered the pairs (xi,yi), ...,{xn,yn) so that Xi < < ... < Then the 

V-statistic for r* is 


iy = \ X! <^(.x^,Xj,Xk,xi)a{yi,yj,yk,yi) 




“ r,,4 I X! 

. l<i<j<k<l<n l<i<j<k<n 


“ r,,4 I X! 

. l<i<j<k<l<n l<i<j<k<n 


^ijkk ^ijjk ^iijk 

2 ^ ^ ^4 

l<i<k<n 

^ijkk '4~ ^ T ^iikk 


E 


l<i<k<n 


Here, the second equality holds since a{xi, Xj,Xki Xi)a{yi, yj,yk, yi) = 0 if any three 
of i, j, k, I are equal. The third equality holds because bijjk = 0 for all i < j < k; 
indeed, Xi < Xj < Xk implies that bijjk corresponds to an inseparable collection of 
points. Note that, in the above equations, we have coefficients of | on bijkkjbujk 
and I on b^kk: these are corrective factors to account for the fact that the number 
of permutations of four elements where exactly two are equal is |S' 4|/2 while the 
number of permutations where exactly two pairs of two are equal is |S'4|/4. Now 
we may continue to rewrite ty as 


~ \ E + E 


^ijkk '4~ 


E 


^iikk 


.l<i<.j<.k<.l<n 


-y, T. 


^ijkl 


l< 2 <_ 7 <fc<n 

^ijkk 


l<i<k<n 


E 


E ^iikl ^iikk 

~9 ^ ^ 


.l<i<j<k<l<n l<i<j<k<n 


l<'i<fc<Z<n 


l<i<k<n 


= E E E E T + E 


^ijkk 


E 


^iikk 


3<k<rL \ fe<Z<n \l^i<j<k 


l<2<fc 


2 41^ 4 

IKiajak l<i</c 


If fc = n then J2k<i<n empty sum which we define to equal 0. For a fixed k < I 
we know already, from Section 3, how to compute J2i<i<j<k efficiently using 
a red-black tree and since bukhbijkk, and b^kk can only correspond to inseparable 












EFFICIENT COMPUTATION OF THE BERGSMA-DASSIOS SIGN COVARIANCE 


11 


or concordant quadruples it is easy to see that 


(A.l) 

(A.2) 

(A.3) 


= 8 ■ {top{k,l)+ bot{k,l)), 


l<i<k 

1 




ftop{k,k)\ /'bot{k,k) 
1 2 )^[ 2 


l<i<j<k 

y^ ^bukk = 4: ■ {top{k,k) + bot{k,k)). 

l<i<k 

Thus we may compute ty by running Algorithm 2 with the following modifications: 
(i) Change line 9 to 

1: for fc = 1 ,n do 


This corresponds to the outer sum of ty. 

(ii) After line 14 add the lines: 

1 : top •<— number of elements > in rbTree 
2 : hot ^ number of elements < yk in rbTree 
3: totalConcordant ^ totalConcordant + ^ ^(* 2 ^ + (^0) 

This accounts for the effect of (A.2) and (A.3). 

(hi) Change line 23 to 

1: totalConcordant *r- totalConcordant + + (*'°*) + ^{top + bot) 

This corresponds to (A.l). 

(iv) Change line 42 to 

1: return ;^(16 • totalConcordant — 8 • total Discordant) 

Finally, note that this Algorithm for computing ty clearly remains 0(n^log(n)). 

Appendix B. Proof of Lemma 1.1 


By permutation invariance, suppose we have relabeled so that xi < X 2 < X 3 < 
X4. We have 3 cases: 


(i) The points in A are inseparable. The fact that &1234 = 0 is an immediate 
consequence of Equation (1.2). 


(ii) The points in A are concordant. In this case we must have that X2 < X3 and 
either max(i/i,y2) < min(y3,y4) or min(?/i,y2) > max(y3,y4). By symmetry 
we need only consider the case when max(?/i,y2) < inin(y3, j/4). By Equation 
( 1 . 2 ) it follows, with some thought, that a(a:,r(i,2,3,4)) = o(j/7r(i,2,3,4)) for all 
permutations ir G S 4 and thus, for any n G S 4 we have n(a::7r(i,2,3,4))a(j/7r(i,2.3,4)) 
a(a^77(i,2,3,4))^ with 


a(®7r(1.2,3,4))^ 


1 if max(a;^(i),a;^( 2 )) < min(a:^( 3 ),a;^( 4 )) or 
min(a:^(i),ai^( 2 )) > max(a:^( 3 ), 0 :^( 4 )) or 
< max(a:^(i), 01 ^( 3 )) < min(a;^( 2 ), 0 :^( 4 )) or 

min(x^(i),a;^(3)) > max(a;^(2),0:^(4)), 

0 otherwise. 


But since xi < X 2 < X 3 < X 4 we have that a(a:: 7 r(i, 2 . 3 , 4 ))a( 2 / 7 r(i, 2 , 3 , 4 )) = 1 if and 
only if { 7 r(l), 7 r( 2 )} e {{1, 2}, {3, 4}} or { 7 r(l), 7 r( 3 )} G {{1, 2}, {3, 4}}. There 
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are exactly 2 ^ = 16 such permutations and thus &1234 = 16. 

(hi) The points in A are discordant. Once again we must have that X2 < X 3 . It 
then follows, by the definition of discordant, that yi ^ and 2/3 ^ j/ 4 . We 
prove an intermediary lemma: 

Lemma B.l. Suppose that (xi, i/i),..., ( 2 : 4 ,2/4) are discordant and xi < X2 < 
X3 < X4. Let 

{x5,y5) = {xi,y2), {xe,ye) = {x2,yi), (2:7,2/7) = (3:^3,2/3), (2^8, ys) = (2:4,2/4), 

so that (2:5,2/5), (2:8,2/8) are simply (2:1,2/1), (2:4,2/4) withyi,y2 switched. 

Then 61234 = 65078- Moreover, the same result is true if we flipped 2 / 3,224 
instead of 2/1,2/2 • 

Proof. First note that, trivially, a( 2 ;,r(i, 2 , 3 , 4 )) = «( 2 : 7 r( 5 , 6 , 7 , 8 )) for ^ ^ 
S4. Let TT be any permutation so that a( 2 ;.n.(i_ 2 , 3 , 4 ))^ = 1- From case (ii) 
we know that we must have { 7 r(l), 7 r( 2 )} £ {{1, 2}, {3,4}} or { 7 r(l), 7 r( 3 )} £ 
{{1,2}, {3,4}}. Suppose that { 7 r(l), 7 r( 2 )} = {1,2}, and let tt' £ 5'4 be the 
permutation where 

7 r'(l) = 7 r( 2 ), 7 r'( 2 ) = 7 r(l), 7 r'( 3 ) = 7 r( 3 ), 7 r'( 4 ) = 7 r( 4 ). 

Then clearly a(2:^(i,2.3,4)) = a(a;,r'(i,2,3,4)) = a(a;,r(5.6.7,8)) = a(2:,7'(5,6.7,8)) but 
a(2/7i-(l,2,3.4)) = a(2/7i-'(5.6,7,8)), a(2/^/(i_2,3,4)) = a(2/7r(5,6,7,8)), 

and thus 

Q(2^7r(l,2,3,4))a(2:7r(l,2,3,4)) + a(2;7r'(1,2,3,4) )o(2:ir'(l,2,3,4) ) 

u(x7r(5,6,7,8) )^(2'7r(5,6,7,8)) 4” n(2:7r'(5,6,7,8) )^(2:7r^(5,6,7,8) ) ■ 

As we may perform a similar procedure to all tt £ 5'4 with a( 2 ;,n.(i 2 , 3 , 4 ))^ = 1 
(changing the choice of tt'), we see that 61234 = 65578 ^ts claimed. 

Finally, pairing tt with tt' given by 

7 r'(l) = 7 r(l), 7 r'( 2 ) = 7 r( 2 ), 7 r'( 3 ) = 7 r( 4 ), 7 r'( 4 ) = 7 r( 3 ) 

shows that this result still holds if we had flipped 2 / 3 , 2/4 instead of 2 / 1 , 2 / 2 - D 

By Lemma B.l, we may assume that xi < 2:2 < 2:3 < X 4 and 2/1 < 2/2 and 
2/3 <2/4- Note that, by the definition of discordant, we must have that 2/2 > 2/3 
and 2/1 < 2/4- From case (ii) we know that there are only 16 permutations tt 
for which a( 2 ; 7 r(i, 2 , 3 , 4 )) 7 ^ 0 and they satisfy 

{7r(l),7r(2)} £ {{1,2}, {3,4}} or {7r(l),^(3)} £ {{1, 2}, {3,4}}. 

If { 7 r(l), 7 r( 2 )} £ {{1,2}, {3,4}} and { 7 r(l), 7 r( 3 )} £ {{1,4}, {2, 3}}, then 
we have a( 2 /,r(i, 2 , 3 , 4 )) = 0. Similarly, a( 2 /,r(i, 2 . 3 , 4 )) = 0 if { 7 r(l), 7 r( 3 )} £ 
{{1,2}, {3,4}} and { 7 r(l), 7 r( 2 )} £ {{1,4}, {2, 3}}. This leaves only 8 per¬ 
mutations TT £ iS '4 for which a( 2 ;,n.(i^ 2 , 3 , 4 ))a( 2 / 7 r(i, 2 , 3 , 4 )) may be non-zero, and 
we check these explicitly: 

a(a:i,2,3,4)0(221,2,3,4) = -1 • 1 = - 1 , 0(212,1,4,3)0(2/2,1,4,3) = -1 • 1 = - 1 , 

o(x 3 ,4,1,2)0(223,4.1,2) = -1 • 1 = - 1 , 0(214,3,2.1)0(2/4.3.2,1) = -1 • 1 = - 1 , 

0(211,3,2.4)0(2/1,3.2.4) = 1 • -1 = - 1 , 0(212,4.1,3)0(2/2.4,1.3) = 1 • -1 = - 1 , 

0(213,1,4.2)0(2/3,1.4,2) = 1 • -1 = - 1 , 0(214,2,3,1)0(2/4.2,3,1) = 1 • -1 = - 1 . 
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We conclude that &1234 = —8 as claimed. 
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